package scu.maqiang.fes;

import java.util.Arrays;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Constants;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES3T101 extends FES {

	public FES3T101(Mesh mesh) {
		super(mesh, 1);
		eleStiff = FES3T101::elementStiff;
		eleMass = FES3T101::elementMass;
		eleSource = FES3T101::elementSource;
	}
	
	public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double[] uu = new double[10];
	    double value = 0;
	    double[][] coord = new double[10][3];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	CG.T10ShapeFunction(coord, xyz, eleVal);
	    	//System.out.println(Arrays.toString(eleVal));
	    	if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
	    		for(int j = 0; j < 10; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		//System.out.println(value);
	    		return value;
	    	}
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	

	public void getValue(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double value = 0;
		int[] ele = mesh.getElement(data.second);
		mesh.getCoordPerElement(ele, coordRegion);
		CG.T6ShapeFunction(coordRegion, xyz, eleVal);
		if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
			data.first = MVO.dot_product(u, ele, eleVal);
			return;
		} else {
		    for(int i = 0, n = mesh.getNt(); i < n; i++) {
		    	ele = mesh.getElement(i);
		    	mesh.getCoordPerElement(ele, coordRegion);
		    	CG.T10ShapeFunction(coordRegion, xyz, eleVal);
		    	//System.out.println(Arrays.toString(eleVal));
		    	if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
		    		for(int j = 0; j < 10; j++) {
		    			eleIdx[j]= ele[j];
		    			value += u[ele[j]] * eleVal[j];
		    		}
		    		data.first = value;
		    		data.second = i;
		    		return;
		    	}
		    }
		    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	   
	}

	/**
	 * ���㵥Ԫ�նȾ���
	 * @param ele
	 * @param xyz
	 * @param kxx
	 * @param tp
	 * @param eleMatrix
	 */
	public static void elementStiff(int[] ele, double[][] xyz, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] N = new double[10];
		double[][] gradN = new double[10][3];
		double[] L = new double[5];
		double area, weight;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 5; i++) {
				L = Constants.TetP3[i].clone();
				area = CG.T10ShapeFunction(xyz, L, N, gradN);
				weight = area * MVO.dot_product(kxx[0], N) * Constants.TetW3[i];
				for (int si = 0; si < 10; si++) {
					for (int sj = 0; sj < 10; sj++) {
						eleMatrix[si][sj] += (gradN[si][0] * gradN[sj][0] + gradN[si][1] * gradN[sj][1] + gradN[si][2] * gradN[sj][2]) * weight;
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for (int m = 0; m < 5; ++m) {
				L = Constants.TetP3[m].clone();
				area = CG.T10ShapeFunction(xyz, L, N, gradN);
				double[][] temp = MVO.matmul(gradN, kxx);
				weight = area * Constants.TetW3[m];
				for (int i = 0; i < 10; i++) {
					for (int j = 0; j < 10; j++) {
						for (int k = 0; k < 3; k++) {
							eleMatrix[i][j] += temp[i][k] * gradN[j][k] * area;
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		
		//System.out.println(MVO.toString(eleMatrix));
		//System.exit(0);
	}

	public static void elementMass(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleMatrix) {
		double[] N = new double[10];
		double[][] gradN = new double[10][3];
		MVO.fill(eleMatrix, 0.0);

		for (int i = 0; i < 5; i++) {
			double[] L = Constants.TetP3[i].clone();
			double area = CG.T10ShapeFunction(xyz, L, N, gradN);
			double weight = area * MVO.dot_product(coef[0], N) * Constants.TetW3[i];
			for (int si = 0; si < 10; si++) {
				for (int sj = 0; sj < 10; sj++) {
					eleMatrix[si][sj] += N[si] * N[sj] * weight;
				}
			}
		}
		if (tp == BVPType.CONSIST_MASS) {
			return;
		} else if (tp == BVPType.LUMP_MASS) {
			double[] sumMatrix = MVO.sumRow(eleMatrix);
			MVO.fill(eleMatrix, 0.0);
			for (int i = 0; i < 10; i++) {
				eleMatrix[i][i] = sumMatrix[i];
			}
		}
	}

	public static void elementSource(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] N = new double[10];
		double[][] gradN = new double[10][3];
		double[] L = null;
		double volume, weight;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 5; i++) {
				L = Constants.TetP3[i].clone();
				volume = CG.T10ShapeFunction(xyz, L, N, gradN);
				weight = volume * MVO.dot_product(f[0], N) * Constants.TetW3[i];
				for (int si = 0; si < 10; si++) {
					eleF[si] += N[si] * weight;
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
//	public void computeGradient(double[] u, double[][] gradU) {
//		int len = u.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		ParamCheck.checkEqual(len, gradU[0].length);
//		ParamCheck.checkEqual(len, gradU[1].length);
//		
//		double[] angle = new double[6];
//		double[][] coord = new double[6][2];
//		double[][] grad = new double[3][2];
//		double[][] eleGradU = new double[6][2];
//		double[] angleCount = new double[this.GetNdof()];
//		Element ele = null;
//		double[] eleU = new double[6];
//		for(int i = 0, n = mesh.getNt(); i < n; i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = u[ele.idx[j]];
//			}
//			CG.gradient3DT4(coord, grad);
//			CG.gradientValue3DT10(grad, eleU, eleGradU);
//			CG.angle2DT6(coord, angle);
//			for(int j = 0; j < 6; j++) {
//				gradU[0][ele.idx[j]] += eleGradU[j][0] * angle[j];
//				gradU[1][ele.idx[j]] += eleGradU[j][1] * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//			//System.out.println(Arrays.toString(MVO.col(eleGradU, 1)));
//		}
//		
//		//System.out.println(MVO.toString(gradU[0]));
//		for(int i = 0, n = mesh.getNv(); i < n; i++) {
//			gradU[0][i] /= angleCount[i];
//			gradU[1][i] /= angleCount[i];
//		}
//	}
//	
//	public void computeSecondDerivative(double[] u, double[][] sdU) {
//		int len = u.length;
//		ParamCheck.checkEqual(len, this.ndof);
//
//		double[] angle = new double[6];
//		double[][] coord = new double[6][2];
//		double[][] grad = new double[3][2];
//		double[] eleSDU = new double[4];
//		double[] angleCount = new double[this.GetNdof()];
//		Element ele = null;
//		double[] eleU = new double[6];
//		for(int i = 0, n = mesh.getNt(); i < n; i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = u[ele.idx[j]];
//			}
//			CG.gradient2DT3(coord, grad);
//			CG.secondDerivative2DT6(grad, eleU, eleSDU);
//			CG.angle2DT6(coord, angle);
//			for(int j = 0; j < 6; j++) {
//				sdU[0][ele.idx[j]] += eleSDU[0] * angle[j];
//				sdU[1][ele.idx[j]] += eleSDU[1] * angle[j];
//				//sdU[2][ele.idx[j]] += eleSDU[2] * angle[j];
//				sdU[3][ele.idx[j]] += eleSDU[3] * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for(int i = 0, n = mesh.getNv(); i < n; i++) {
//			sdU[0][i] /= angleCount[i];
//			sdU[1][i] /= angleCount[i];
//			//sdU[2][i] /= angleCount[i];
//			sdU[3][i] /= angleCount[i];
//		}
//		
//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[0][ele.idx[j]];
//			}
//			sdU[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[3][ele.idx[j]];
//			}
//			sdU[3][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[3][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[3][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//		}
//		
//		for(int i = 0; i < mesh.getNv(); i++) {
//			sdU[2][i] = sdU[1][i];
//		}
//	}
//
	public void computeFirstAndSecondDerivative(double[] u, double[][] fdU, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		double[][] coord = new double[10][3];
		double[][] grad = new double[4][3];
		double[][] eleGradU = new double[10][3];
		double[] eleSDU = new double[10];
		double[] volumeCount = new double[this.GetNdof()];
		double[] eleU = new double[10];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 10; j++) {
				eleU[j] = u[ele[j]];
			}
			double volume = CG.gradient3DT4(coord, grad);
			CG.gradientValue3DT10(grad, eleU, eleGradU);
			CG.secondDerivative3DT10(grad, eleU, eleSDU);
			for(int j = 0; j < 10; j++) {
				fdU[0][ele[j]] += eleGradU[j][0] * volume;
				fdU[1][ele[j]] += eleGradU[j][1] * volume;
				fdU[2][ele[j]] += eleGradU[j][2] * volume;
				sdU[0][ele[j]] += eleSDU[0] * volume;
				sdU[1][ele[j]] += eleSDU[1] * volume;
				sdU[2][ele[j]] += eleSDU[2] * volume;
				sdU[3][ele[j]] += eleSDU[3] * volume;
				sdU[4][ele[j]] += eleSDU[4] * volume;
				sdU[5][ele[j]] += eleSDU[5] * volume;
				volumeCount[ele[j]] += volume;
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			fdU[0][i] /= volumeCount[i];
			fdU[1][i] /= volumeCount[i];
			fdU[2][i] /= volumeCount[i];
			sdU[0][i] /= volumeCount[i];
			sdU[1][i] /= volumeCount[i];
			sdU[2][i] /= volumeCount[i];
			sdU[3][i] /= volumeCount[i];
			sdU[4][i] /= volumeCount[i];
			sdU[5][i] /= volumeCount[i];
		}
		
		for(int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			for(int k = 0; k < 6; k++) {
				for(int j = 0; j < 10; j++) {
					eleU[j] = sdU[k][ele[j]];
				}
				sdU[k][ele[4]] = 0.5 * (eleU[0] + eleU[1]);
				sdU[k][ele[5]] = 0.5 * (eleU[1] + eleU[2]);
				sdU[k][ele[6]] = 0.5 * (eleU[2] + eleU[3]);
				sdU[k][ele[7]] = 0.5 * (eleU[0] + eleU[3]);
				sdU[k][ele[8]] = 0.5 * (eleU[1] + eleU[3]);
				sdU[k][ele[9]] = 0.5 * (eleU[2] + eleU[3]);				
			}
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[6]] = 0.5 * (eleU[2] + eleU[3]);
//			sdU[1][ele.idx[7]] = 0.5 * (eleU[0] + eleU[3]);
//			sdU[1][ele.idx[8]] = 0.5 * (eleU[1] + eleU[3]);
//			sdU[1][ele.idx[9]] = 0.5 * (eleU[2] + eleU[3]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[6]] = 0.5 * (eleU[2] + eleU[3]);
//			sdU[2][ele.idx[7]] = 0.5 * (eleU[0] + eleU[3]);
//			sdU[2][ele.idx[8]] = 0.5 * (eleU[1] + eleU[3]);
//			sdU[2][ele.idx[9]] = 0.5 * (eleU[2] + eleU[3]);
		}
	}

}
